Exponentiated Weibull distribution (exponweib)#
The exponentiated Weibull distribution (SciPy: scipy.stats.exponweib) is a flexible continuous distribution on \([0,\infty)\) (in its standard form) obtained by raising the Weibull CDF to a positive power.
It is popular in reliability and survival analysis because—with two shape parameters—it can represent a wide range of hazard-rate behaviors (increasing, decreasing, unimodal, and bathtub-like).
Learning goals#
Write down the PDF/CDF/PPF and connect them to the Weibull distribution.
Understand how the parameters change shape and hazard behavior.
Derive moments using a clean change-of-variables + series expansion.
Implement inverse-CDF sampling using only NumPy.
Use SciPy (
exponweib) for evaluation, simulation, and MLE fitting.
import numpy as np
import plotly
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
import scipy
from scipy import integrate, optimize, special
from scipy.stats import exponweib as exponweib_dist
from scipy.stats import weibull_min, chi2, kstest
# Plotly rendering (CKC convention)
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
# Reproducibility
rng = np.random.default_rng(42)
np.set_printoptions(precision=6, suppress=True)
# Record versions for reproducibility (useful when numerical details matter)
VERSIONS = {"numpy": np.__version__, "scipy": scipy.__version__, "plotly": plotly.__version__}
VERSIONS
{'numpy': '1.26.2', 'scipy': '1.15.0', 'plotly': '6.5.2'}
1) Title & Classification#
Name:
exponweib(Exponentiated Weibull; SciPy:scipy.stats.exponweib)Type: Continuous
Support (standard form): \(x \in [0,\infty)\)
Parameter space (standard form): shape parameters \(a>0\), \(c>0\)
SciPy location/scale:
loc \in \mathbb{R},scale > 0with\[X = \text{loc} + \text{scale}\,Y, \qquad Y \sim \mathrm{ExponWeib}(a,c).\]
Unless stated otherwise, we work with the standard form (loc=0, scale=1).
Notation note: many texts write the parameters as \((\alpha, k, \lambda)\), where \(\alpha\) is the exponentiation parameter, \(k\) is the Weibull shape, and \(\lambda\) is the scale. SciPy uses \((a, c)\) for the two shape parameters and
scalefor \(\lambda\).
2) Intuition & Motivation#
What it models#
Start with a Weibull random variable \(W\) (shape \(c\), scale \(\lambda\)):
The exponentiated Weibull distribution raises this CDF to a power \(a>0\):
So it is literally the Weibull CDF, exponentiated.
A useful interpretation when \(a\) is an integer (\(a=m\in\mathbb{N}\)):
If \(W_1,\dots,W_m\) are i.i.d. Weibull\((c,\lambda)\), then $\(\max\{W_1,\dots,W_m\} \sim \mathrm{ExponWeib}(a=m, c, \lambda).\)$
This “maximum-of-\(m\) Weibulls” story helps build intuition: larger \(a\) pushes mass to the right (stochastically larger lifetimes).
Typical real-world use cases#
Reliability engineering: modeling lifetimes with non-monotone hazard (e.g. bathtub-shaped failure rates).
Survival analysis: flexible parametric baseline hazard.
Hydrology / environmental extremes: positive-valued quantities with skew and flexible tails.
Manufacturing / materials: time-to-failure with early-life and wear-out effects.
Relations to other distributions#
\(a=1\) gives the standard Weibull distribution (
scipy.stats.weibull_min).\(c=1\) gives the exponentiated exponential distribution.
\(a=1\) and \(c=1\) gives the exponential distribution.
Location/scale transform recovers different units and origins: \(X=\text{loc}+\text{scale}\,Y\).
3) Formal Definition#
We present the standardized form first (loc=0, scale=1).
CDF#
For \(x \ge 0\), \(a>0\), \(c>0\):
PDF#
For \(x>0\):
Quantile function (PPF)#
For \(q\in(0,1)\):
This inverse-CDF is the key to fast sampling.
Location / scale#
With loc and scale>0, define \(z = (x-\text{loc})/\text{scale}\). Then for \(x\ge \text{loc}\):
Hazard function (survival analysis)#
The hazard rate is
For the exponentiated Weibull, the hazard can be increasing, decreasing, unimodal, or bathtub-shaped depending on \((a,c)\).
def exponweib_cdf(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""CDF of the exponentiated Weibull (NumPy implementation).
Matches SciPy's `scipy.stats.exponweib.cdf` for the same (a, c, loc, scale).
"""
if a <= 0 or c <= 0 or scale <= 0:
raise ValueError("Require a>0, c>0, scale>0")
x = np.asarray(x, dtype=float)
z = (x - loc) / scale
out = np.zeros_like(z, dtype=float)
mask = z >= 0
zm = z[mask]
# exm1c = 1 - exp(-zm**c), computed stably for small zm**c
exm1c = -np.expm1(-(zm**c))
out[mask] = exm1c**a
return out
def exponweib_logpdf(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""Log-PDF of the exponentiated Weibull (stable for small/large x)."""
if a <= 0 or c <= 0 or scale <= 0:
raise ValueError("Require a>0, c>0, scale>0")
x = np.asarray(x, dtype=float)
z = (x - loc) / scale
out = np.full_like(z, -np.inf, dtype=float)
mask = z > 0
zm = z[mask]
neg_zc = -(zm**c)
exm1c = -np.expm1(neg_zc) # = 1 - exp(-zm**c) in (0,1)
logp = (
np.log(a)
+ np.log(c)
- np.log(scale)
+ (a - 1.0) * np.log(exm1c)
+ neg_zc
+ (c - 1.0) * np.log(zm)
)
out[mask] = logp
return out
def exponweib_pdf(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""PDF computed from the log-PDF."""
return np.exp(exponweib_logpdf(x, a, c, loc=loc, scale=scale))
def exponweib_ppf(q: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""Percent-point function (inverse CDF)."""
if a <= 0 or c <= 0 or scale <= 0:
raise ValueError("Require a>0, c>0, scale>0")
q = np.asarray(q, dtype=float)
if np.any((q < 0) | (q > 1)):
raise ValueError("q must be in [0,1]")
out = np.empty_like(q, dtype=float)
out[q == 0] = loc
out[q == 1] = np.inf
mask = (q > 0) & (q < 1)
qm = q[mask]
out[mask] = loc + scale * (-np.log1p(-(qm ** (1.0 / a)))) ** (1.0 / c)
return out
def exponweib_rvs_numpy(
a: float,
c: float,
*,
loc: float = 0.0,
scale: float = 1.0,
size=1,
rng: np.random.Generator | None = None,
) -> np.ndarray:
"""Draw samples from ExponWeib(a,c) using **only NumPy**.
Inverse transform sampling:
U ~ Uniform(0,1)
X = loc + scale * (-log(1 - U^{1/a}))^{1/c}
"""
if a <= 0 or c <= 0 or scale <= 0:
raise ValueError("Require a>0, c>0, scale>0")
rng = np.random.default_rng() if rng is None else rng
u = rng.uniform(size=size)
# Avoid log(0) at the endpoints.
eps = np.finfo(float).eps
u = np.clip(u, eps, 1.0 - eps)
return loc + scale * (-np.log1p(-(u ** (1.0 / a)))) ** (1.0 / c)
def exponweib_sf(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""Survival function S(x)=1-CDF(x), computed stably via log1p."""
cdf = exponweib_cdf(x, a, c, loc=loc, scale=scale)
return np.exp(np.log1p(-cdf))
def exponweib_hazard(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
"""Hazard h(x)=f(x)/S(x), computed in log-space for stability."""
x = np.asarray(x, dtype=float)
logf = exponweib_logpdf(x, a, c, loc=loc, scale=scale)
sf = exponweib_sf(x, a, c, loc=loc, scale=scale)
out = np.zeros_like(x, dtype=float)
mask = sf > 0
out[mask] = np.exp(logf[mask] - np.log(sf[mask]))
out[~mask] = np.inf
return out
# Sanity check: our formulas match SciPy (prefer logpdf to avoid overflow near 0).
a, c, loc, scale = 1.7, 1.3, 0.4, 2.2
x = np.logspace(-5, 2, 40) * scale + loc
qs = np.linspace(0.01, 0.99, 9)
assert np.allclose(exponweib_logpdf(x, a, c, loc=loc, scale=scale), exponweib_dist.logpdf(x, a, c, loc=loc, scale=scale))
assert np.allclose(exponweib_cdf(x, a, c, loc=loc, scale=scale), exponweib_dist.cdf(x, a, c, loc=loc, scale=scale))
assert np.allclose(exponweib_ppf(qs, a, c, loc=loc, scale=scale), exponweib_dist.ppf(qs, a, c, loc=loc, scale=scale))
# Check that inverse-CDF sampling produces the right distribution (quick KS test).
x_samp = exponweib_rvs_numpy(a, c, loc=loc, scale=scale, size=3_000, rng=rng)
D, p = kstest(x_samp, lambda t: exponweib_dist.cdf(t, a, c, loc=loc, scale=scale))
print(f"KS statistic={D:.4f}, p-value={p:.4f}")
KS statistic=0.0092, p-value=0.9576
4) Moments & Properties#
Raw moments#
Assume the standard form with loc=0 and scale=\lambda.
Using the substitution \(y = (x/\lambda)^c\) (so \(x=\lambda y^{1/c}\)), the exponentiated Weibull density simplifies nicely and we obtain a series expression for the raw moments:
For integer \(a=m\in\mathbb{N}\), the binomial series terminates at \(j=m-1\).
For non-integer \(a\), the series is infinite but often converges quickly.
From the raw moments we get:
mean: \(\mu = \mathbb{E}[X]\)
variance: \(\sigma^2 = \mathbb{E}[X^2] - \mu^2\)
skewness and kurtosis from the central moments \(\mu_3,\mu_4\).
MGF / characteristic function#
There is no simple elementary closed form for the MGF in general. Existence depends on the tail:
If \(c>1\), the tail is super-exponential and \(M_X(t)=\mathbb{E}[e^{tX}]\) exists for all real \(t\).
If \(c=1\) (exponentiated exponential tail), \(M_X(t)\) exists for \(t < 1/\lambda\).
If \(0<c<1\), \(M_X(t)\) diverges for any \(t>0\) (but the Laplace transform exists for \(t\le 0\)).
The characteristic function \(\varphi_X(t)=\mathbb{E}[e^{itX}]\) always exists and can be computed numerically.
Entropy#
The differential entropy is
SciPy provides exponweib.entropy(...) (numerical integration). A Monte Carlo estimate using -mean(logpdf(samples)) is also useful for sanity checks.
def exponweib_raw_moment_series(
r: float,
a: float,
c: float,
*,
scale: float = 1.0,
max_terms: int = 50_000,
tol: float = 1e-12,
) -> float:
"""Compute E[X^r] via the binomial/Gamma series (loc=0).
Uses the identity:
E[X^r] = scale^r * a * Gamma(1 + r/c) * sum_{j>=0} (-1)^j * binom(a-1, j) / (j+1)^{1 + r/c}
Notes
-----
- Requires a>0, c>0, scale>0.
- For integer a, the sum is finite (j=0,...,a-1).
- For non-integer a, the sum is infinite; we truncate by `tol`.
"""
if a <= 0 or c <= 0 or scale <= 0:
raise ValueError("Require a>0, c>0, scale>0")
if r <= -c:
raise ValueError("Moment diverges for r <= -c in the standard form")
p = 1.0 + r / c
prefactor = (scale**r) * a * special.gamma(1.0 + r / c)
s = 0.0
for j in range(max_terms):
coeff = special.binom(a - 1.0, j)
term = ((-1.0) ** j) * coeff / (j + 1.0) ** p
s_new = s + term
# Stop once the marginal contribution becomes negligible.
if j > 0 and abs(term) < tol * max(1.0, abs(s_new)):
s = s_new
break
s = s_new
return float(prefactor * s)
def exponweib_mvsk_from_raw_moments(a: float, c: float, *, scale: float = 1.0):
"""Return mean, variance, skewness, excess kurtosis via raw moments."""
m1 = exponweib_raw_moment_series(1.0, a, c, scale=scale)
m2 = exponweib_raw_moment_series(2.0, a, c, scale=scale)
m3 = exponweib_raw_moment_series(3.0, a, c, scale=scale)
m4 = exponweib_raw_moment_series(4.0, a, c, scale=scale)
mean = m1
var = m2 - m1**2
mu3 = m3 - 3 * m2 * mean + 2 * mean**3
mu4 = m4 - 4 * m3 * mean + 6 * m2 * mean**2 - 3 * mean**4
skew = mu3 / (var ** 1.5)
kurt_excess = mu4 / (var**2) - 3.0
return mean, var, skew, kurt_excess
# Compare the series-based moments to SciPy's numerical stats.
a, c, scale = 1.6, 1.2, 2.0
mean_s, var_s, skew_s, kurt_s = exponweib_mvsk_from_raw_moments(a, c, scale=scale)
mean_sp, var_sp, skew_sp, kurt_sp = exponweib_dist.stats(a, c, scale=scale, moments="mvsk")
print("Series moments:")
print(" mean=", mean_s)
print(" var =", var_s)
print(" skew=", skew_s)
print(" kurt(excess)=", kurt_s)
print("\nSciPy stats (numerical):")
print(" mean=", float(mean_sp))
print(" var =", float(var_sp))
print(" skew=", float(skew_sp))
print(" kurt(excess)=", float(kurt_sp))
Series moments:
mean= 2.4294714656111736
var = 2.7333895473411074
skew= 1.3022259281402193
kurt(excess)= 2.474111034724598
SciPy stats (numerical):
mean= 2.429471462970749
var = 2.7333895579577634
skew= 1.30222591791057
kurt(excess)= 2.4741110129320223
# Entropy: SciPy (numerical integration) vs Monte Carlo.
a, c, scale = 1.6, 1.2, 2.0
rv = exponweib_dist(a, c, scale=scale)
H_scipy = rv.entropy()
x = rv.rvs(size=60_000, random_state=rng)
H_mc = -np.mean(rv.logpdf(x))
print(f"Entropy (SciPy) = {H_scipy:.6f}")
print(f"Entropy (Monte Carlo)= {H_mc:.6f}")
Entropy (SciPy) = 1.756397
Entropy (Monte Carlo)= 1.757041
# Characteristic function (CF) example via Monte Carlo
a, c, scale = 1.2, 0.9, 1.5
rv = exponweib_dist(a, c, scale=scale)
x = rv.rvs(size=120_000, random_state=rng)
def cf_mc(t: float) -> complex:
return complex(np.mean(np.exp(1j * t * x)))
for t in [0.5, 1.0, 2.0, 4.0]:
print(f"t={t:>3}: phi(t)≈ {cf_mc(t)}")
t=0.5: phi(t)≈ (0.5638539622591024+0.489223667427389j)
t=1.0: phi(t)≈ (0.2521834895613472+0.4322972648554105j)
t=2.0: phi(t)≈ (0.07024447912956258+0.2711293646917929j)
t=4.0: phi(t)≈ (0.01435532670293679+0.14062825155057773j)
5) Parameter Interpretation#
Think of the parameters as controlling (i) the underlying Weibull shape and (ii) how strongly we “exponentiate” its CDF.
Shape parameter \(c\) (Weibull shape)#
\(c<1\): heavy mass near 0, decreasing hazard (infant mortality)
\(c=1\): constant hazard (exponential-like tail)
\(c>1\): increasing hazard (wear-out)
Exponentiation parameter \(a\)#
Because \(F(x) = F_W(x)^a\):
\(a>1\) makes the distribution stochastically larger (shifts probability mass to the right).
\(0<a<1\) makes it stochastically smaller (more mass near 0).
For integer \(a=m\), \(X\) is the maximum of \(m\) i.i.d. Weibull lifetimes.
Scale and location#
scalerescales the x-axis (units of time/length/etc.).locshifts the support to start atloc.
Below we visualize how the PDF, CDF, and hazard change with parameters.
def plot_pdf_grid(param_sets, *, title: str):
fig = go.Figure()
for (a, c, scale) in param_sets:
rv = exponweib_dist(a, c, scale=scale)
x_max = rv.ppf(0.995)
x = np.linspace(1e-6, x_max, 500)
fig.add_trace(go.Scatter(x=x, y=rv.pdf(x), mode="lines", name=f"a={a}, c={c}, scale={scale}"))
fig.update_layout(title=title, xaxis_title="x", yaxis_title="pdf")
fig.show()
def plot_cdf_grid(param_sets, *, title: str):
fig = go.Figure()
for (a, c, scale) in param_sets:
rv = exponweib_dist(a, c, scale=scale)
x_max = rv.ppf(0.995)
x = np.linspace(0.0, x_max, 500)
fig.add_trace(go.Scatter(x=x, y=rv.cdf(x), mode="lines", name=f"a={a}, c={c}, scale={scale}"))
fig.update_layout(title=title, xaxis_title="x", yaxis_title="cdf")
fig.show()
def plot_hazard_grid(param_sets, *, title: str):
fig = go.Figure()
for (a, c, scale) in param_sets:
rv = exponweib_dist(a, c, scale=scale)
x_max = rv.ppf(0.995)
x = np.linspace(1e-6, x_max, 600)
h = rv.pdf(x) / rv.sf(x)
fig.add_trace(go.Scatter(x=x, y=h, mode="lines", name=f"a={a}, c={c}, scale={scale}"))
fig.update_layout(title=title, xaxis_title="x", yaxis_title="hazard h(x)")
fig.show()
# Vary 'a' at fixed c
plot_pdf_grid([(0.6, 1.5, 1.0), (1.0, 1.5, 1.0), (2.0, 1.5, 1.0)], title="PDF: effect of a (fixed c=1.5, scale=1)")
plot_cdf_grid([(0.6, 1.5, 1.0), (1.0, 1.5, 1.0), (2.0, 1.5, 1.0)], title="CDF: effect of a (fixed c=1.5, scale=1)")
# Vary 'c' at fixed a
plot_pdf_grid([(1.5, 0.7, 1.0), (1.5, 1.0, 1.0), (1.5, 2.0, 1.0)], title="PDF: effect of c (fixed a=1.5, scale=1)")
plot_hazard_grid([(1.2, 0.7, 1.0), (1.2, 1.0, 1.0), (1.2, 2.0, 1.0)], title="Hazard: effect of c (fixed a=1.2, scale=1)")
# A few parameter combinations that show different hazard shapes
plot_hazard_grid(
[(0.6, 0.8, 1.0), (1.0, 1.0, 1.0), (2.0, 0.9, 1.0), (2.0, 2.0, 1.0)],
title="Hazard: different shapes across (a, c)",
)
6) Derivations#
We sketch derivations for the expectation/variance and the likelihood.
6.1 Expectation (raw moments)#
Work with loc=0, scale=\lambda and consider the \(r\)-th raw moment:
With
use the substitution \(y=(x/\lambda)^c\) (so \(x=\lambda y^{1/c}\) and \(dx = (\lambda/c) y^{1/c-1}\,dy\)). The Jacobian cancels the Weibull power term and the integral becomes
Now expand using the binomial series (valid for \(0<e^{-y}<1\)):
Swap sum and integral (justified under standard conditions) to get
Finally, recognize a Gamma integral:
which yields the moment formula used earlier.
6.2 Variance#
Compute \(\mathbb{E}[X]\) and \(\mathbb{E}[X^2]\) from the raw-moment expression and combine:
6.3 Likelihood (iid sample)#
Given iid data \(x_1,\dots,x_n\) from the standard model with scale \(\lambda\) and shapes \((a,c)\), the likelihood is
The log-likelihood (often optimized numerically) is
SciPy’s fit routine maximizes this (with loc/scale included) via numerical optimization.
def exponweib_loglik(x: np.ndarray, a: float, c: float, *, loc: float = 0.0, scale: float = 1.0) -> float:
"""Total log-likelihood using our NumPy logpdf."""
return float(np.sum(exponweib_logpdf(x, a, c, loc=loc, scale=scale)))
def fit_exponweib_mle_via_minimize(x: np.ndarray, *, loc_fixed: float = 0.0):
"""Simple MLE demo using SciPy optimize on transformed parameters.
We optimize over (log a, log c, log scale) with loc fixed.
"""
x = np.asarray(x, dtype=float)
def nll(theta):
log_a, log_c, log_scale = theta
a = np.exp(log_a)
c = np.exp(log_c)
scale = np.exp(log_scale)
return -exponweib_loglik(x, a, c, loc=loc_fixed, scale=scale)
# crude initial guess from a Weibull fit
c0, loc0, scale0 = weibull_min.fit(x, floc=loc_fixed)
theta0 = np.log([1.0, c0, scale0])
res = optimize.minimize(nll, theta0, method="Nelder-Mead")
a_hat, c_hat, scale_hat = np.exp(res.x)
return a_hat, c_hat, loc_fixed, scale_hat, res
# Quick demo on synthetic data
true = dict(a=1.8, c=1.1, loc=0.0, scale=2.5)
x = exponweib_dist.rvs(true["a"], true["c"], loc=true["loc"], scale=true["scale"], size=900, random_state=rng)
(a_hat, c_hat, loc_hat, scale_hat, res) = fit_exponweib_mle_via_minimize(x, loc_fixed=0.0)
print("True params:", true)
print("Minimize MLE:", {"a": a_hat, "c": c_hat, "loc": loc_hat, "scale": scale_hat})
print("Converged:", res.success)
True params: {'a': 1.8, 'c': 1.1, 'loc': 0.0, 'scale': 2.5}
Minimize MLE: {'a': 1.6943268223977832, 'c': 1.1384941757974947, 'loc': 0.0, 'scale': 2.556153488212525}
Converged: True
7) Sampling & Simulation#
NumPy-only algorithm (inverse transform)#
Because the CDF has a closed-form inverse, sampling is straightforward.
Let \(U\sim\mathrm{Uniform}(0,1)\). Set
Then \(X \sim \mathrm{ExponWeib}(a,c,\text{loc},\text{scale})\).
Implementation detail: for numerical stability,
compute \(1-\exp(-t)\) as
-expm1(-t)compute \(\log(1-u)\) as
log1p(-u)clip \(U\) away from exactly 0 or 1 to avoid
log(0).
# Sampling with the NumPy-only sampler
a, c, scale = 1.6, 1.2, 2.0
x = exponweib_rvs_numpy(a, c, scale=scale, size=60_000, rng=rng)
# Compare empirical moments to theory
mean_emp = x.mean()
var_emp = x.var()
mean_theory, var_theory = exponweib_dist.stats(a, c, scale=scale, moments="mv")
print(f"Empirical mean={mean_emp:.4f}, variance={var_emp:.4f}")
print(f"Theory mean={float(mean_theory):.4f}, variance={float(var_theory):.4f}")
Empirical mean=2.4353, variance=2.7418
Theory mean=2.4295, variance=2.7334
8) Visualization#
We visualize:
the PDF for different parameter settings
the CDF (as a sanity check for probability mass)
Monte Carlo samples vs the theoretical PDF
# PDF / CDF and Monte Carlo comparison
a, c, scale = 1.6, 1.2, 2.0
rv = exponweib_dist(a, c, scale=scale)
x_max = rv.ppf(0.995)
xs = np.linspace(1e-6, x_max, 600)
fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=rv.pdf(xs), mode="lines", name="pdf"))
fig.update_layout(title="Exponentiated Weibull PDF", xaxis_title="x", yaxis_title="pdf")
fig.show()
fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=rv.cdf(xs), mode="lines", name="cdf"))
fig.update_layout(title="Exponentiated Weibull CDF", xaxis_title="x", yaxis_title="cdf")
fig.show()
# Monte Carlo check
samples = rv.rvs(size=40_000, random_state=rng)
fig = go.Figure()
fig.add_trace(go.Histogram(x=samples, nbinsx=60, histnorm="probability density", name="samples", opacity=0.7))
fig.add_trace(go.Scatter(x=xs, y=rv.pdf(xs), mode="lines", name="theoretical pdf"))
fig.update_layout(title="Samples vs theoretical PDF", xaxis_title="x", yaxis_title="density")
fig.show()
9) SciPy Integration#
SciPy exposes this distribution as scipy.stats.exponweib with methods:
pdf,logpdf,cdf,sf,ppfrvsfor random variatesstats/moment/entropyfitfor maximum-likelihood estimation
A common workflow:
Pick a parametric family (here,
exponweib).Fit it to data.
Check fit visually (histogram, QQ plot) and via likelihood-based criteria.
Use the fitted model for inference or simulation.
# SciPy usage: pdf/cdf/rvs/fit
# 1) Create a frozen distribution
rv = exponweib_dist(1.6, 1.2, loc=0.0, scale=2.0)
xs = np.linspace(0, rv.ppf(0.99), 6)
print("x grid:", xs)
print("pdf:", rv.pdf(xs))
print("cdf:", rv.cdf(xs))
# 2) Generate data
x = rv.rvs(size=2_000, random_state=rng)
# 3) Fit parameters (MLE)
# For strictly-positive lifetime data it's often sensible to fix loc=0.
a_hat, c_hat, loc_hat, scale_hat = exponweib_dist.fit(x, floc=0.0)
print("\nFitted params (floc=0):")
print({"a": a_hat, "c": c_hat, "loc": loc_hat, "scale": scale_hat})
# 4) Compare to a Weibull fit (nested model a=1)
a0, c0, loc0, scale0 = exponweib_dist.fit(x, fa=1.0, floc=0.0)
print("\nWeibull-as-exponweib (fa=1) fit:")
print({"a": a0, "c": c0, "loc": loc0, "scale": scale0})
x grid: [0. 1.548109 3.096218 4.644327 6.192436 7.740545]
pdf: [0. 0.295522 0.171121 0.069922 0.024508 0.00785 ]
cdf: [0. 0.351985 0.721415 0.899531 0.967212 0.99 ]
Fitted params (floc=0):
{'a': 1.388222664766818, 'c': 1.2931256583561437, 'loc': 0.0, 'scale': 2.2045665712074545}
Weibull-as-exponweib (fa=1) fit:
{'a': 1.0, 'c': 1.5427975513751417, 'loc': 0.0, 'scale': 2.7025618683212835}
10) Statistical Use Cases#
10.1 Hypothesis testing (nested models)#
Because the Weibull distribution is the special case \(a=1\), you can test
\(H_0: a=1\) (Weibull) vs
\(H_1: a \ne 1\) (exponentiated Weibull)
using a likelihood ratio test (LRT):
10.2 Bayesian modeling#
Bayesian modeling treats \((a,c,\lambda)\) as random and combines a prior with the likelihood:
There is no conjugacy here, but generic MCMC (e.g. Metropolis-Hastings) works well.
10.3 Generative modeling#
Once fitted (frequentist or Bayesian), the distribution can be used to:
simulate lifetimes for stress testing
generate synthetic positive-valued data with realistic skew/tails
build parametric simulators inside larger pipelines
# 10.1 Likelihood ratio test: Weibull (a=1) vs Exponentiated Weibull
# Simulated data under the alternative
true = dict(a=1.8, c=1.1, loc=0.0, scale=2.5)
x = exponweib_dist.rvs(true["a"], true["c"], loc=true["loc"], scale=true["scale"], size=900, random_state=rng)
# Fit alternative (a free) with loc fixed to 0
(a1, c1, loc1, scale1) = exponweib_dist.fit(x, floc=0.0)
ll1 = float(np.sum(exponweib_dist.logpdf(x, a1, c1, loc=0.0, scale=scale1)))
# Fit null (a=1) with loc fixed to 0
(a0, c0, loc0, scale0) = exponweib_dist.fit(x, fa=1.0, floc=0.0)
ll0 = float(np.sum(exponweib_dist.logpdf(x, a0, c0, loc=0.0, scale=scale0)))
lr_stat = 2.0 * (ll1 - ll0)
p_value = chi2.sf(lr_stat, df=1)
print("Alt fit:", {"a": a1, "c": c1, "scale": scale1})
print("Null fit:", {"a": a0, "c": c0, "scale": scale0})
print(f"LR statistic={lr_stat:.3f}, p-value={p_value:.4g}")
Alt fit: {'a': 2.2710907010507073, 'c': 0.9742363049697826, 'scale': 2.082356652940917}
Null fit: {'a': 1.0, 'c': 1.4663239585093428, 'scale': 3.7688940378192104}
LR statistic=21.856, p-value=2.94e-06
# 10.2 Bayesian modeling (toy example): random-walk Metropolis on log-parameters
x = exponweib_dist.rvs(1.8, 1.1, scale=2.5, size=500, random_state=rng)
# Prior: independent normals on log-parameters (broad, weakly informative)
prior_mu = np.array([0.0, 0.0, 0.0])
prior_sigma = np.array([1.5, 1.5, 1.5])
def log_prior(theta: np.ndarray) -> float:
d = theta.size
return float(
-0.5 * np.sum(((theta - prior_mu) / prior_sigma) ** 2)
- np.sum(np.log(prior_sigma))
- 0.5 * d * np.log(2 * np.pi)
)
def log_lik(theta: np.ndarray) -> float:
log_a, log_c, log_scale = theta
a, c, scale = np.exp([log_a, log_c, log_scale])
return float(np.sum(exponweib_dist.logpdf(x, a, c, loc=0.0, scale=scale)))
def log_post(theta: np.ndarray) -> float:
return log_lik(theta) + log_prior(theta)
n_steps = 6_000
burn = 1_500
step_scale = np.array([0.08, 0.08, 0.08])
# Initialize at the MLE (good starting point)
a_mle, c_mle, loc_mle, scale_mle = exponweib_dist.fit(x, floc=0.0)
cur = np.log([a_mle, c_mle, scale_mle])
cur_lp = log_post(cur)
chain = np.empty((n_steps, 3), dtype=float)
accepted = 0
for t in range(n_steps):
prop = cur + step_scale * rng.normal(size=3)
prop_lp = log_post(prop)
if np.log(rng.uniform()) < (prop_lp - cur_lp):
cur, cur_lp = prop, prop_lp
accepted += 1
chain[t] = cur
acc_rate = accepted / n_steps
print(f"Acceptance rate: {acc_rate:.3f}")
post = chain[burn:]
a_s, c_s, scale_s = np.exp(post.T)
summary = {
"a_mean": float(a_s.mean()),
"c_mean": float(c_s.mean()),
"scale_mean": float(scale_s.mean()),
"a_95%": tuple(np.quantile(a_s, [0.025, 0.975])),
"c_95%": tuple(np.quantile(c_s, [0.025, 0.975])),
"scale_95%": tuple(np.quantile(scale_s, [0.025, 0.975])),
}
summary
Acceptance rate: 0.178
{'a_mean': 1.4587089908602129,
'c_mean': 1.2533464400252319,
'scale_mean': 3.084652557279298,
'a_95%': (1.0627486616378665, 2.2450850386814643),
'c_95%': (1.0083703414202325, 1.4947448812498432),
'scale_95%': (2.1649431209050984, 3.8530214219000483)}
# Simple posterior trace plots
fig = go.Figure()
fig.add_trace(go.Scatter(y=a_s, mode="lines", name="a"))
fig.update_layout(title="Posterior trace: a", xaxis_title="iteration", yaxis_title="a")
fig.show()
fig = go.Figure()
fig.add_trace(go.Scatter(y=c_s, mode="lines", name="c"))
fig.update_layout(title="Posterior trace: c", xaxis_title="iteration", yaxis_title="c")
fig.show()
fig = go.Figure()
fig.add_trace(go.Scatter(y=scale_s, mode="lines", name="scale"))
fig.update_layout(title="Posterior trace: scale", xaxis_title="iteration", yaxis_title="scale")
fig.show()
# 10.3 Generative modeling: posterior predictive vs fitted MLE
# Use the same x from the Bayesian section
# MLE predictive
rv_mle = exponweib_dist(a_mle, c_mle, loc=0.0, scale=scale_mle)
x_mle = rv_mle.rvs(size=15_000, random_state=rng)
# Posterior predictive (draw parameters, then draw a sample)
idx = rng.integers(0, len(a_s), size=400)
pp_samples = []
for i in idx:
rv_i = exponweib_dist(a_s[i], c_s[i], loc=0.0, scale=scale_s[i])
pp_samples.append(rv_i.rvs(size=40, random_state=rng))
pp_samples = np.concatenate(pp_samples)
fig = go.Figure()
fig.add_trace(go.Histogram(x=x, nbinsx=60, histnorm="probability density", name="observed", opacity=0.6))
fig.add_trace(go.Histogram(x=x_mle, nbinsx=60, histnorm="probability density", name="MLE predictive", opacity=0.5))
fig.add_trace(go.Histogram(x=pp_samples, nbinsx=60, histnorm="probability density", name="Posterior predictive", opacity=0.5))
fig.update_layout(title="Observed vs predictive distributions", barmode="overlay", xaxis_title="x", yaxis_title="density")
fig.show()
11) Pitfalls#
Invalid parameters: require \(a>0\), \(c>0\),
scale>0.Behavior at \(x\approx 0\): depending on \(c\) and \(a\), the density may go to 0 or blow up; use
logpdffor stability.CDF cancellation: for very small \(x\), \(1-\exp(-x^c)\) suffers cancellation; use
expm1(-expm1(-t)).Inverse-CDF sampling: avoid exact 0/1 uniforms; use
clipandlog1p.Fitting: unconstrained
loccan drift negative even for strictly-positive data; fixfloc=0when appropriate.Moment calculations: the binomial series can require many terms for extreme parameters; prefer SciPy numerical moments or Monte Carlo when in doubt.
MGF existence: for \(c\le 1\), the MGF may not exist for positive \(t\); use the characteristic function or Laplace transform instead.
12) Summary#
exponweibis a continuous distribution on \([0,\infty)\) (standard form) with shape parameters \(a>0\) and \(c>0\).It is defined by exponentiating a Weibull CDF: \(F(x) = [1-\exp(-x^c)]^a\) (with optional
loc/scale).It generalizes Weibull (
a=1) and can represent diverse hazard shapes important in reliability/survival analysis.Raw moments admit a useful Gamma + binomial series representation; SciPy provides numerical
statsandentropy.Sampling is easy via inverse transform and can be implemented with NumPy only.
References#
SciPy documentation:
scipy.stats.exponweib.Wikipedia: “Exponentiated Weibull distribution”.